Disentangling multi-level systems: averaging, correlations and memory 

Jeroen Wouters, Valerio Lucarini 

Meteorologisches Institut, 
University of Hamburg, 
KlimaCampus 
Grindelberg 7, 20144 Hamburg, Germany 

October 28, 2011 

Abstract 

We consider two weakly coupled systems and adopt a perturbative approach based on the Ruelle response 
theory to study their interaction. We propose a systematic way to parametrize the effect of the coupling as a 
function of only the variables of a system of interest. Our focus is on describing the impacts of the coupling 
on the long-term statistics rather than on the finite-time behaviour. By direct calculation, we find that, at first 
order, the coupling can be surrogated by adding a deterministic perturbation to the autonomous dynamics of 
the system of interest. At second order, there are additionally two separate and very different contributions. 
One is a term taking into account the second order contributions of the ffuctuations in the coupling, which 
can be parametrized as a stochastic forcing with given spectral properties. The other one is a memory term, 
coupling the system of interest to its previous history, through the correlations of the second system. If these 
correlations are known, this effect can be implemented as a perturbation with memory on the single system. In 
order to treat this case, we present an extension to Ruelle's response theory able to deal with integral operators. 
We discuss our results in the context of other methods previously proposed to disentangle the dynamics of two 
coupled systems. We emphasize that our results do not rely on assuming a time scale separation, and, if such a 
separation exist, can be used equally well to study the statistics of the slow as well as that of the fast variables. 
By recursively applying the technique proposed here, we can treat the general case of multilevel systems. 

1 Introduction 

The dynamics of complex systems can often be partitioned in such a way to highlight a multi-level structure, 
where each level consists of a subset of variables that interact strongly with each other and are weakly coupled 
with variables belonging to other levels of the system. Often, the dynamics of each level takes place in distinct 
spatial and temporal scales, so that it is possible to hierarchically arrange the variables from the most to the 
least relevant, in view of the properties of interest. In practical applications, the observation or the simulation 
at high resolution of processes occurring on small temporal and spatial scales is very challenging. It is crucial 
to devise strategies for accounting, at least in an approximate way, the impact of the fast, almost inaccessible 
variables, onto the slow variables, which are more usually the main object of investigation. In the language of 
dynamical systems theory, the overarching goal is to define rigorously or empirically an effective dynamics for the 
slow manifold of the system able to take into account the influence of the, enslaved, fast manifold. Otherwise, 
the presence of such very different scales make these systems hard to simulate, so that they are often referred 
to as stiff problems (H], and hard to observe in practical terms, because no overall optimal spatial and temporal 
resolution for the collection of data can be deflned, with the result that some variables are poorly described, 
whereas other variables have redundant observations. 

In these terms, the climate system provides an especially interesting case study. In the study of geophysical 
fluid dynamics, one often encounters a wide range of scales. The spatial scales involved in atmospheric dynamics 
can range from lO^m down to 10~^m. Temporal scales can vary from lO^^s to 10~®s [161 130j . Here, one has to 
take into account coupling between different subsystems, which feature dynamics taking place in vastly different 
ranges, as clarifled when considering separately the atmosphere, the oceans, the biosphere, and the criosphere 
[211 128] . We still do not have a standard formulation of a climate model, because, depending on the class of 
problems we want to tackle, we need to resort to different approximations and even different formulations of 
the relevant dynamics. The state-of-the-art climate models used for the IPCC projections (2^ cannot be used 
to study the last million years of Earth's history, not only because they are too computationally expensive, but 
because they do not take into the right considerations the processes which dominate the Earth dynamics on such 



long time scales. On the other hand, Earth Models of Intermediate complexity which are tailored for 

such problems, are definitely not competitive if our aim is to describe the variability of the present climate and 
the short-term climate change. The problem of providing a satisfactory description of the dynamical processes 
that cannot be explicitly be represented has long been addressed at various levels of rigour by the climate 
community [201 1311 1111 1101 [S] 1281 122) . The solutions that are being used in practical applications today go 
by the name of parametrizations. Processes that are not explicitly described are approximated via functions 
(deterministic or, more recently, stochastic) of the variables of the model, basically along the lines of classical 
closure theories in fluid dynamics for the Reynolds' stress tensor. Given the approximate and empirical nature 
of the parametrizations, a problem that often arises is their robustness with respect to changes in the considered 
modelling set-up, such as the case when additional variables are included in the model or different large scale 
conditions are considered. Retuning the parametrization is definitely a delicate aspect of climate modelling 
and a major hurdle towards the goal of obtaining flexible climate models able to obtain the so-called seamless 
simulation of its statistical properties and seamless (probabilistic) prediction of its state. 

In this paper, we tackle this class of problems from a rather general point of view. We will consider a two-level 
system, consisting of two weakly coupled systems X and Y . We focus our attention on the X system. We want 
to derive ab initio a suitable parametrization of its coupling with the Y system. 

For this purpose, various reduction techniques have been proposed in the literature. A classic approach is 
the averaging method O [T^ , which assumes that there is a vast time-scale separation between the two systems 
X and Y . Such method allows to derive a renormalised dynamics for the X system whose trajectories converge 
on finite time scales to those of the original dynamics. An interesting use of the averaging method for deriving a 
simplified dynamics has been proposed in [1 , where a sort of mean field ansatz is taken. In statical mechanics, 
projector operator techniques which allow to formally eliminate certain variables have found a great deal of 
success, as in the case of the Mori-Zwanzig equation [33) . 

In the present manuscript, we take a different route, based on the response theory developed by Ruelle [^l27l 
114) . Such a theory provides explicit formulas for computing how the long-time averages of Axiom A dynamical 
systems (described by the so called SRB measure [32]) change as a result of small perturbations to the flow. This 
approach is especially useful for studying the impact of changes in the internal parameters of a system or of small 
modulations to the external forcing [241 [5l [2l 1151 119[ . According to the theory, the changes in the observables 
can be expressed basically as a power series in orders of the strength of the forcing, so that one should assume, 
in our case, that the coupling between the X and Y systems - which acts as perturbation to the dynamics of X 
- should be small. We derive explicit expressions for the perturbations in the statistics of the X system due to 
the coupling at the first and second order. We explain how these can be parametrized in terms of a surrogate 
deterministic forcing, multiplicative noise, and memory term, so that we construct an approximated projected 
dynamics on the X variables. Our approach does not make use of any large time scale separation between the 
X and the Y variables, and does not even assume that the X variables are slower than the Y variables. 

This article is structured as follows. In Section[2]we give a short review of Ruelle's response theory at higher 
orders. In Section [3] we then use this theory to derive the response of the system X to the coupling with the 
system Y . We derive the first and second order terms and present an expression for a surrogate dynamics of X 
alone. The correction to the unperturbed dynamics will include a deterministic forcing, a stochastic forcing, and 
a forcing given by a term which takes into account the history of the system. We present also a diagrammatic 
representation of the most relevant contributions. In Section |4] we briefly discuss our results in comparison to 
classic variables reduction methods such as the averaging method [31 E] and the projection operator method 
[33) . In Section [5] we present our conclusions and perspectives for future work. In Appendix |A] we describe how 
to extend the Ruelle response formula to forcings describing memory effects, which we use to treat the memory 
terms resulting at the second order of perturbation. 

2 Response theory 

Ruelle [26) recently derived explicit formulas for describing the smooth dependence of the SRB measure of 
Axiom A dynamical systems to small perturbations of the flow. Such response theory boils down to a Kubo- 
like perturbative expression connecting the terms describing the linear and nonlinear response of the system 
as expectation values of observables on the unperturbed measure. At order n of nonlinearity, such expectation 
values can be written as n— uple convolution of a causal Green function with the time-delayed perturbative 
fields, so that at every order Kramers-Kronig relations can be written for the Fourier transform of the Green 
function, the so-called susceptibility '14'. In order to clarify the goal of response theory and set our notation, we 
present a simple derivation of the linear and higher-order nonlinear response formulas obtained by introducing 
the so-called Dyson interaction picture for the operators. Given a continuous time dynamical system Xt £ M 
with dynamics 

xt = F{xt) (1) 



2 



and a corresponding flow function /' relating Xt to the initial conditions xq, Xt — f*{xo), we assume that an 
invariant measure po exists ^32,, describing long-term averaged quantities: 

po{A) = hm U^dT [ l{dx)A{r{x)) . (2) 

We want to know how such averages, given by the so called Sinai- Ruelle-Bowen (SRB) invariant measure po, 
are affected by small changes in the dynamics of the system. Let's assume that the dynamics of the system is 
altered by adding a time independent perturbation ^^(a::): 

Xt = F{xt) + *(a;t) , (3) 

with /' the perturbed flow function, so that we end up with a modified SRB measure p: 



p{A) = ^lim - dr j l{dx)A{rix)) . (4) 

If the strength of the perturbation is small, the corresponding change in the SRB measure of the system will 
also in general be small. The occurence of resonances has been studied in |27l SI |6] , but this is not expected to 
occur in system of physical interest (see [7] and references therein). Hence, we can expand the perturbed SRB 
measure in powers of the perturbation ^I*. In order to do so, we first expand the evolved observable A{ft{x)). 
This amounts to constructing a Dyson series, where the perturbation to the dynamical system plays the role of 
the interaction Hamiltonian. Using the chain rule, the time-dependence of an observable A at time t is described 
by 

^^it^ = CA{f\xo)) = {Co + Cr)A{f\xo)) , (5) 

where 

CoA{x)^F{x)-V^A{x), £i = ■ V.A(x) , C = Co + Ci. (6) 

By defining an interaction observable as Ai{t, xq) — Ilo{~T)A{f^{xo)) — exp(— £01") exp{jCT)A{xo) with IIo{t) = 
exp{£oT), we can make an expansion of A in orders of the perturbation 4'. This observable Aj obeys the 
differential equation 

dAi(T,xo) r t \ A f ^ /'7\ 
— = jCi(t)Ai(t, xo) , (7) 

where 

LiiT) = Uo{-r)Cino{T). (8) 
Note that the operator Ci is of order ^. A formal solution to Eq. [7]is given by Duhamel's formula: 

Ai{t,xo) = Ai{0,xo)+ [ dsiCi{si)Ai{si,xo) . (9) 
Jo 

By iterating this solution and making use of the fact that Ai{t, xo) = exp{— CoT)A{f'^ (xq)) , we get an expression 
for the evolved observable A{f^{xo)) in terms of A{0,xo): 



Aifixo)) =no{T)A{xo) + / dsiHoir - s^)£iUo{si)A{xo) 

Jo 

+ [ dsi [ ' dS2llo(T - Si)£ino(si - S2)£ino(s2)^(xo) + 0(*^) . 

Jo Jo 

Such an expression for the expectation values under perturbed dynamics A{f'^{xo)) in terms of the unperturbed 
dynamics I1o{t)A{xo) — A{f'^{xo)) allows us to expand the SRB measure p around po- The perturbed SRB 
measure p can be written as 

p{A) = ^lirn^ i ^ dr y l{dx) (no{T)A{x) + ^ ds,Ilo{T - si)JZiIlo{si)A{x) + 0{<f'')^ 

= po(A) + 5p(^)(^) + 0(*^), (10) 



3 



where the first term in this sum gives the unperturbed measure po, while the second term is the linear response 
to '5. We can rewrite it as 



dsi J dr J l{dx)Il„{T - si)CiTlo{si)A{x) 

1 fT [T-^i r 
lim — / ds\ I du I l{dx)Il(,{u)C\Iio{si)A{x) 



5pW(^)= hm ' '''' 







poo poo p 

/ ds^poiC^Uo{s^)A) = ds^ po{dx)^ix) ■ \/^A(nx))) 
Jo Jo J 

dsiG%{s,). (11) 



This shows that the linear response can be written in terms of averages taken in the unperturbed system. 
Carrying on the calculation, we obtain the following expression for the higher order terms: 

poo poo /"^l 

piA) = po{A)+ dsipoiCiUo{si)A) + dsi ds2po(£ino(si)£ino(s2)^) + 0(*^) 
Jo Jo Jo 

= po{A) + Sp^'^A) + 5p^^\A) + 0(*^) , (12) 

where 

poc f^l 

Sp^^\A)= dsi ds2po{Cino{si)Cino{s2)A) 
Jo Jo 

poo poo p 

= / dsi / ds2 po{dx)<5>ix)-\/,^{rix))-\/fs^^,^A{r+'^x))) 
Jo Jo J 

i poc 

dsi / ds2G^'^(si, S2) . (13) 
^0 

3 Two-level systems 

We will now make use of the response formulae derived in Section [2] to study how the statistical properties of a 
dynamical system which is composed, in the unperturbed case, by two independent subsystems, change when a 
weak coupling is introduced. Let's assume that the variables x of the dynamical system can be separated into 
two groups of variables X and Y 



Y ^ 

so that the unperturbed dynamical system can be considered as given by the uncoupled dynamics: 

X^Fx{X), 

Y^Fy{Y). (14) 
A weak coupling is then introduced between the two subsystems as follows: 

X ^ Fx{X) + <fxiX,Y), 

Y^Fy{Y) + 9y{X,Y). (15) 
The unperturbed and perturbed evolution operators Co and £1 are now given by 

Co = Fx{X)-Vx + Fy{Y)-Vy (16) 

and 

£1 = *x(X, F) ■ Vx + *y (X, Y) ■ Vy 

= *^-V, (17) 

where = *y) and = (Vx Vy). 
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3.1 First order 

We derive from Eqs. [6] and |12| that the first order response to £i is given by 

roo 

= / dsipoi£ino{si)A) (18) 
Jo 

roo r 

= j dsi J po{dx)'i>''{x)\/{Ao f^){x), (19) 

where po is the SRB measure corresponding to the uncoupled dynamics of Eq. |14| The response is an integral 
over time of the response function GJ;^^ (r) describing the delayed impact of a forcing in the past over a time 
span of r: 

p oo 

5«p(^) = / dsiG«(si) 
Jo 

G%{si) = J po{dx)<f^{x)-\/{Aor^){x). (20) 

By exploiting the fact that the unperturbed dynamics of the X and Y subsystems are uncoupled, this becomes 

G%{r) = I Poidx)^^{x) ■ VA (^IfyJ) ) , (21) 

where fx and fy represent the unperturbed flows of the X and Y subsystems, respectively. We are only 
interested in the evolution of observables on the X system so we consider observables A that depend on X 
variables only: 



) = A(X) . 



(22) 



Therefore, the response function G^'^ can be written as 



G%{r) ^ j PO *x(X,y) . VxAi fxiX)) . (23) 

Since the dynamics of the unperturbed system is decoupled, its SRB measure po can be written as a product of 
the SRB measure po^x{dX) and po,Y{dY). We have: 



G%(r) = j po,x(dX){*x(X, Y))p„.^ ■ VxAifxiX)) . (24) 

Such a first order perturbation can be represented in a diagrammatic way as depicted in Figure [l] where the 
arrow indicates the direction of the interaction and the time variable refers to the presence of a time-delayed, 
causal effect on the system of interest of the external perturbation. 

Equation [24] can be used to derive a parametrization of the first order effect of coupling the X and Y 
subsystems on the A{X) observables, expressed as a function of the X variables only. In other terms, we want 
to derive an expression for a perturbation vector field Mi such that if we consider the dynamical system 

^=Fx(X) + Mi(X), (25) 



where Fx{X) is the same as in Eq. 14 up to the first order the expectation values of the A observables over 
the SRB invariant measure pi of the dynamical system given in Eq. [25] and over the invariant measure px,Y are 
identical. The invariant measure po{dX) for the unperturbed Mi — system is, obviously, identical to partial 
measure po.x{dX) introduced in Eq. 23 Therefore, we must impose that the first order effect in Eq. 25 due to 
All is the same as the first order correction due to the coupling between the X and the Y variables introduced 
in Eq. [15] We discover that by choosing: 

Mi(X) :={*x(X,r))p„ ,, (26) 
and plugging this expression in Eq. |11| we obtain that up to the first order in 'i' and VA: 

Sp^^l{A)^5i'^p(A), (27) 

such that 

Px.y{A) = po,x{A) + Si'^piA) + 0(vl>^) = pi (A) + 0(*=') = po(A) + 5p« (A) + 0(*2) . (28) 

Concluding, the effect on X of its coupling with Y can be described at first order with a mean-field theory 
constructed by taking the F-ensemble average of the coupling, where the averaging is taken on the uncoupled 
Y dynamics. 
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3.2 Second order 

From Eq. |l2] we can derive that at second order the response is a double integral of a response function G^'^ 
by plugging the expressions given in Eqs. |16| and [TtI into Eq. [13] 



S. 



POG POO POO POO 

= / dsi / ds2 p{CiIlo{si)CiIlois2)A) = / dsi / ds2G^_^^^^(si, S2) 
Jo JO JO JO 

= 4%.aP(^) + 4%,i,P(^) = / ds, dS2G^%,^,,(si,S2)+G^%_^,,(si,S2), (29) 



('2') (1') 

where we separate such a second order contribution into the sum of two terms, G\ ^, ^ ^ and G\ ^, ^ ^, defined 
as 

Gi'*,*,a(si,S2) = j po(dx)*x(X,y)- Vxno(si)(*x(X, Y)-Vx){A{r{X)), 
G^'**.(,(si.S2) = / Po{dx)<i>Y{X,Y)-VYllo{si){^x{X,Y)-Vx){A{r{X)). (30) 



The Green function G)j \j ^ a('''i' '''2) responsible for second contribution 5^ \j, ^ can be depicted in a diagrammatic 
mode as in Figure |2] This term describes the second order effect of the perturbation of the Y dynamics on the 
statistical properties of the X system, and does not take into account the modifications on the statistical proper- 
ties of Y due to its coiling with X. Instead, Figure 3 provides a graphical representation of the G^'^ ^ i,(Ti, T2), 

responsible for the second order contribution ^ j,. In this case, the second order effect comes as a result of the 
modification to the first order effect shown in Figure [l]due to the impact of the X dynamics on the Y systems: 
note that in this case, as opposed to Figures [T][2j one arrow points from the X to the Y system. 

Our goal is now to provide a second order parametrization of the coupling in terms of the X variables alone, 
so that the invariant measure p2 of the resulting system is such that VA: 

px.y(A) = p2(yl) + 0(*'). (31) 

In the following, we show that can be accomplished by the following reduced system: 

X(t) = Fx(X(t)) + Afi(X(t)) + ^*;,,(X(t))a,(t)+ / drhir, X{t - r)) 

= Fx{X{t)) + Kh{X{t)) + M2{X{t)) + K'h{X(t)) , (32) 

where the first term on the right hand side given the unperturbed dynamics, the second term is the same given 
in the first order approximation presented in Eq. |25[ the third term is given by a combination of stochastic 
forcings, and the fourth term describes a memory effect. For simplicity, we indicate all of the perturbative terms 
as M-j(X(t)) wi th an abuse of notation. We will now provide the explicit expression for the third and fourth 
terms in Eq. l32|by matching the corrections up to 0($^) to the unperturbed dynamics of the previous with what 
given in Eqs. |20| and Eq. |29[ Obviously, since Eq. [32] has a stochastic component, its invariant measure f>2 will 
be accordingly defined. Please note that in the following, when we refer to first or second order perturbations 
to the measure po due to the Mj terms in Eq. |32[ the expectation value over the realizations of the stochastic 
processes are implicitly considered. Note that in |17) we have explained how to treat the impact of stochastic 

ical systems along the lines of the Ruelle response theory. 



perturbations on deterministic dynan 



The first order response in Eq. 



32 



can be written as (5^^|^p(y4), where 



roc n nCD 

5^mAA) = / dsi / poMdX) / dsi G^\jsi) 
Jo J Jo 

= dsi J po,x{dX)MkiX)-VxA{fjl{X)). (33) 



The full second order response can be written as ^ 5\l ^,j.p(A), where 



(34) 



POO POO P POO POO 

'5m1,m,p(^) = / dsi I ds2 / po,x{dX) j dsi / ds2 G^^^j^-u.mA'^i, ^2) 
Jo Jo J Jo Jo 

roc pec p 

= ds^ ds2 J po.x{dX)Nh{X) - yxMAr^iX)) - V ^s,^^^A{r^^+'-{X)) . 
Combining Eqs. |28| and [3T] we impose that 

4"p(^) + 4%,.p(^) + 4%,.p(^) = E <p(^) + E 4'!.M,p(^) + o(*') ■ (35) 
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Our procedure will lead to simplifying Eq. [35] into the following: 

Si'^piA) + 5i%_^p{A) + = S'ilpiA) + S[%p{A) + S^^l,,^p{A) + ^Iai.pC^) , (36) 



where the ter ms missing on the right hand side in Eq. 36 with respect to Eq. 35 are of order 0(^'^) or higher. 

From Eq. we have that S'^^ p{A) = 
observe that \j, ^p{A) can be written as: 



From Eq. 28 we have that 5^^^ p{A) = We now consider the next order of perturbation. We first 



4%,aP(^) = / dsi ds2 / po(dx)((*x(X,y) - M{X)) ■ Vx)no(si)((*x(X,y) - Mix)) ■ Vx)iA{f'X) 
Jo Jo J 

poo roo r 

+ dsi ds2 po{dx){M{X)-Vx)Ilo{si){(A4{X))-Vx){A{f^X) 
Jo Jo J 

= I po{dx){<S>'x{X,Y) ■Vx)no{s,){^'x{X,Y) ■Vx)A{nX) + S^l^,,^p{A) 

= 4%,^p(A) + <5g)_,,^p(A), (37) 

where 'ii'x{X,Y) — 'iix{X,Y) — M{X) gives the fluctuations of ^' around its Y averaged value M. We now 
construct the term M2{X{t)) so that, i1 

M2^'(^) corresponds exactly to 5^'^ fP{A) introduced in Eq. 



construct the term M2{X(t)) so that, its first order impact (5^^^p(A) vanishes, whereas its second contribution 

which we refer to as the fluctuation term. 



37 



37 



we flrst consider the case where the 



Fluctuation term 

To analyze the S^^^ fp{A) = J J ds-ids2G^^^ ^ f{s\,S2) term in Eq 
coupling fluctuations 'il'x are separable in the X and Y variable: 

Vx{X,Y) = -i!'x,i{X)'^'x,2iX) ■ 
In this case the response function G\\^ ^ ^ can new written as: 

G%,^A'^^'^) = I Po,x(dX)3(si)*:^,i(X)ax*:^,i(/xU^))9;n(;,)(Ao/i-)(/ii(X)), (38) 

where 

g{si) = {*k2(y)*x,2(P (y)))po,. ■ (39) 
We now impose that the second order contributions 5^'^ fP{A) and S^^^^ j^j-^p{A) are identical, so that: 

'^A,*,*,/('*i' ^a) = G'^%^2,M2(*i' ■'2) 

= J po,xidX)M2{T - s, - S2,X)dxM2iT - s^, nX))df.,^x)iAir'+'^X)) (40) 

for some time-dependent perturbation A'h- This can be achieved by choosing the perturbation M2 to be a 
multiplicative noise term described as follows: 

M2(t,X) = *:^,i(X)a(t), 
{a{ti)a{t2)) = g{t2 - ti) , 

{a{t))=0, (41) 

where the average in the second equality is an average over realizations of the noise. Taking the average of the 
model response over the realizations of the noise gives 

{Gam.,M2 (*!' ^2)) = G%,^j{si,S2) => 5i|,M,p(A) = 4%,/p(^) • (42) 

These results can be easily generalized also to the case where the coupling term 'i/'x is not separable. In this 
case, 'I'^ can be decomposed in product of Schauder bases of function in X and Y: 

*:y(x,y) = ^*:^,i,,(x)vi>:^,2,,(y). 



The response function is now given by 



Ga 



%j{si,S2)^ I Po,xidX) J2 9,Asi)^'xMX)dx<l''x^MfxiX))d^^^,{Aorx-){rx^{X)) 

i,j,k,l 

gjAsi)^{^'x,2AY)^'x,2AfYiY)). (43) 



7 



The response can in this case be modelled with a linear combination of multiplicative noise 

M2(X) = ^vI>;,,(X)a,(t), (44) 

{a,{h)ai{t2)}=g,,,(t2-ti), (45) 
{a,it))^0. (46) 

Note also that, since {uj(t)) — Vj, we have that 5^j^j^p{A) = 5^^' J^^^p(y4) — 5^^^ j^j-^p{A) = 0, where we 
have already anticipated that Ms is not written as a stochastic term. At this stage, we have shown that 

Memory term 

In order to complete the construction of the parametrization leading to the surrogate dynamics given in Eq. |36| 
we need to show that 5^'^ bPi^) ~ ^M^Pi'^)^ ^^'^ that the two terms (5^J' m^p{^) a^nd (5^^^ MsPC^) 0{'^^) or 



(2) 

higher. The G\\^ ^ j, term in Eq. 30 can be seen as a modification to the first order response on the X system 
related to the coupling on Y due to the fact that the Y system is modified by the interaction with the X system. 
After integrating over the Y variables, we can express 5^^^^, f,p{A) as: 

pCG poo r 

/ ds, / ds2 po.x{dX){^YAX,Y)dY,.<l>xAr{X),r{Y))},,_^df.,^x),AA(r+''{X)). 
Jo Jo J 

(47) 

This can be expressed as: 

roc POO f 

ds^ ds2 / Po,x(dX)/i,(s2,X)a/=2(x),,(^(r^+""(X)), (48) 
Jo Jo J 

where we introduce the following integral kernel: 

hAs,x) = {*y,.(x,y)ay,,*x,,(r(x),r(y)))p„,^ . (49) 

This term describes a memory effect: the evolution of the Y variables between times T — si — S2 and T ~ si 
keeps track (at first order) of the previous state of the X system. Therefore, when we consider the impact of 
the Y variables on the X system, at second order we get a contribution coming from the history of X itself. 
Obviously, the faster the decorrelation time of the Y system, the weaker the influence of this term. 

A last step needed at this stage is to provide a response formula for integral forcings, since the formulas 
derived in Section[2]are only valid for perturbations without memory. In Appendix|A]we derive such a generalized 
expression, which includes as a special case that which has been discussed in Section [2] The linear response to 
a perturbation which can be written as follows: 

roc 

X{t) = Fx{X{t)) ^ X{t) = Fx{X{t)) + drh{t-r,X{t^r)) (50) 

Jo 

results to be V observables A: 

pCG poo p 

S^'A{A)= ds^ ds2 po,x{dX)h,{s2,X)dfs,t,x)AM.n^"AX)). (51) 
Jo Jo J 



This indeed implies that by choosing Afs = drh{T, X{t — r)) in Eq. 32 we obtain bPi-^) ~ ^^M^Pi-^)- 
is important to note that, by construction, even if 5'^^^p{A) is a first order perturbation of the expectation value 



of ^, it is 0{^^). Therefore, all terms such as 5^^^ ^^p(A), 5^^j_^ j^^^p^A) presented in Eq. 35 are 0(5"^) or higher. 



so that Eq. |36]is indeed correct. Analogously to the case of the noise discussed in the previous section, the term 
hj (s, X) can be estimated from a trajectory of the system and then used to provide the suitable parametrization. 
This concludes our derivation and leads to the result anticipated in Eq. |32| 

3.3 Independent coupling case 

An especially interesting result can be obtained when considering the not uncommon situation where the coupling 
is independent of the variable it is affecting, i.e. 'i!x{Y) depends only on Y and ^y{X) depends only on X. 
In this case, the reduced dynamical system able to mimic the statistics of the full system for the observables 
depending on X variables only up to the second order on the coupling strength can be written as: 

pOQ 

X(s)=Fx{X(s)) + Mi+a{s)+ dr h{s ~ r, X' {t - r)) . (52) 

Jo 
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The first order term Mi is now a constant drift term and the fluctuation term is an additive noise term with 
correlations 

Hh)a{t2)) = - Mi)(*x(/'^-'^ {¥)) - Kh)),, ,. . (53) 

The integral kernel of the memory term simplifies to 

h{s,X) = *y(X)(9y*x(r(V)))po,v ■ (54) 

where the role of the impact of the memory effect of the Y variables on the memory of the X system is especially 
clear. 

4 Comparison with other variables reduction methods 

The problem of variables reduction is a classic topic in dynamical systems theory and statistical mechanics. We 
will briefly discuss the relation between these methods and what presented here. A well-developed approach is 
the so called averaging method jSjilSj. This method deals with coupled systems of the form 

X = f(X,Y), (55) 

Y=-^g(X,Y). (56) 

In the limit where the time-scale difference between the two goes to infinity, by taking e — >■ oo, the trajectories 
of the original system converge on finite time scales to the trajectories of an averaged version of the X system, 
namely 

X' = py|x(/(X',y)) , (57) 

where Y is average over the measure Py\x, which is the SRB measure corresponding to the dynamics of Y for 
fixed X: 

Y^gx{Y), (58) 
gx{Y) ^ g(X,Y) . (59) 

The fluctuations around this average trajectory can be approximated by an Ornstein-Uhlenbeck process. One 
should note that while we have discussed here long term averages in the form of SRB measures, the averaging 
method describes convergence on finite time scales. Another point of difference is the weak coupling assumption 
we have made here. This assumption is not made in the averaging method, where the limit of an infinitely fast 
Y process is taken instead. Such a limit basically takes care of getting rid of the higher order effects discussed 
in this paper. 

In [l] the averaging method is used to derive a simplified reduced model. The X dependent measure pY\x is 
first replaced by the measure Py\x' for a fixed value X* for X. To improve the statistics of the reduced model 
a linear term in X is added: 

X' = pY^x {f{X', Y)) + {X- X*)5pYix' {f(X', Y)) , (60) 

where SpY\x* is the response of the measure Py\x' to a change in the value of the fixed X. Such a feedback 
mechanism, where a change in X affects the X dynamics through the presence of Y, is also present in the memory 
effect we find here. However, as the limit of infinite time-scale differences is taken in the averaging method, there 
the feedback happens instantly from the point of view of the slow process and no memory is present. 

Projector operator techniques have found much use in the statistical mechanics for deriving Fokker-Planck 
equations. A celebrated result in the Mori-Zwanzig equation, which consists of an integro-differential equation, 
containing an integral term describing the memory effect due to projected out variables (33,. Consider the 
setting of a dynamical system with Liouville operator L where we have a projector operator P projecting out 
irrelevant parameters. Every observable A can the be then be written as a two-vector containing Ai = PA and 
A2 = QA = (1 — P)A. The Liouville equation for this vector is then given by 



dt 1^2/ ^21 L22) \A2 



(61) 



where Ln = PLP, L12 = PLQ, L21 = QLP and L22 = QLQ. By using Duhamel's formula for the equation 
for A2 a solution for A2 can be found where the term with L21A1 is seen as a perturbation to the free evolution 
of A2 under L22 • By plugging this solution back into the differential equation for Ai , the following equation of 
motion for Ai is obtained 

d /■* 

—Ai{t) = LiiAi{t) + L12 J dsexp{L22{t-s))L2iAi{s) + Li2exp{L22t)A2{0). (62) 
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Again there is a feedback mechanism through the presence of the irrelevant dynamics. It has to be noted that 
the Mori-Zwanzig equation is not an approximation, but an equivalent formulation of the original problem. 
This approach is concerned with describing the instantaneous evolution of probability measures, describing the 
evolution of the current state, rather than the long term averages, as discussed here. Our strategy may provide 
a way to compute, at least approximately the projector operators that in the Mori-Zwanzig theory are defined 
only formally. 

5 Conclusions 

In this paper we have adopted the Ruelle response theory to study the effect of a weak coupling between two 
dynamical systems. We have derived explicit formulas to compute the changes in the expectation value of a 
general observable of one of the two systems (we call this system X, and refer to the other system as Y) up to 
the second order in the coupling strength. We have also presented a simple diagrammatic representation of such 
terms, which is suggestive of the possibility of reframing the Ruelle theory in graphical terms along the lines of 
the Feynman representation of quantum electrodynamics. 

These results allow us to define a surrogate perturbed dynamics for the system X such up to second order 
its statistical properties are the same as those of the coupled system when a general observable of X alone is 
considered. In other terms, we have attempted to derive a parametrization of the coupling. We have found that 
the invariant measure for the coupled system restricted to the X system can be approximated by the invariant 
measure of a reduced system with three forcing terms Mi: 

j^X{t) = Fx{X{t)) + Mi(X(f)) + M2{X{t)) + M3(X(t)), (63) 

where Fx{X{t)) gives the unperturbed flow. Mi is a deterministic correction, which fully describes the change 
in the invariant measure up to first order, M2 is a stochastic forcing expressed as a sum of multiplicative noise 
terms with non-trivial correlation properties, and M3 is an integral expression which describes a memory term. 
The explicit expressions for the Mj terms are given in Eq. |32[ The Mi term basically describes a mean field 
theory, where the effective X dynamics is modified by adding the effect of the Y averaged dynamics. The M2 
term, instead, takes care of representing the effect of the fluctuations of the Y systems. Finally, the Mg term 
describes the non-markovian effect due to the X system impacting the state of the Y state at a given time in 
the past, which in turn perturbs the X system at a later time. Note that, if the coupling is independent of the 
variable it is affecting, so that the perturbation on X (Y) depends only on Y {X) variables, the parametrization 
becomes especially attractive, as Mi reduces to a constant, Ai2 to an additive noise term and M3 has a greatly 
simplified integral kernel. 

It is important to note that in our derivation no assumptions have been made on the dynamics of X and Y 
systems in terms of differing time scales and/or of decorrelation time, our only hypothesis being the assumption 
that the coupling is weak in some sense: there is - a priori - symmetry between the X and Y systems, we just 
select for which one we want to derive an effective dynamics. Actually, the characteristics of the unperturbed 
dynamics of the two systems are expected to play an important role in the relative importance of the fiuctuation 
and memory terms. In a system that is highly mixing, correlations in time disappear fast and the fiuctuation 
term will act as white noise. The memory term on the other hand will quickly drop off to zero in this case. If 
the Y system decorrelates slowly, there will be a more pronounced memory effect. 

It is also important to note that our method does not aim at describing within a certain accuracy and 
for a certain time horizon the individual orbits of the system of interest, but rather focuses on providing an 
approximate method for reconstructing its statistical properties. Therefore, the method discussed in this paper 
has points of contact and differences with variable reduction techniques previously proposed in the literature 
such as the averaging method, which assumes a very large time scale separation but no limit on the strength 
of the coupling, and the projection operator method which instead is more formal and has less straightforward 
practical applicability. 

Obviously, if the system consists of more than two levels it can be reduced to only one by applying repeatedly 
to the approach described in this paper, so that our method has a good degree of generality when variables 
reduction is sought. 

These findings will be of particular interest in the scientific fields dealing with multilevel systems and in 
particular, in the context of climate modelling. First of all we have derived an explicit strategy to substitute 
unresolved processes that is robust to changes in parameters or scales. Furthermore are the parameters (such 
as the mean coupling, the noise correlations and the memory kernel) explicitly formulated in terms of statistical 
properties only of the unresolved processes. By comparing the fiuctuation and memory terms we see that 
different strategies will be applicable whether the unresolved process is fast or slow. Hence, this is relevant when 
modelling, e.g., the impact of the atmosphere on the oceans versus the other way around. 

As future research we will be implementing the proposed modelling strategies in low-dimensional systems, 
such as the Lorenz '96 system [13] . Hopefully this can lead to applications on more realistic models. Of particular 
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interest will be how well this derivation based on weak couplings will apply to couplings in the climate system 
and whether the interplay between the fluctuations and memory can be observed. 
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A Response to perturbations with memory 

As in [18], we derive the response of a discrete dynamical system 

xt+i = f{xt) ■ 

Consider a perturbation x that has a memory of the variable at I time steps in the past: 

xt+i = f{it) +xi{xt-i)- 
We can expand evolution of x around the evolution of x. By induction, one can show that 

k 

xt+k = /"{it) + J2 Df''-'{f{xt))xi{xt-i+,.i) + O(x') , (64) 
j=i 

or on an observable A that 

k 

A{x^+k) = Aifii,)) + J2DiAo f''-'){f{xt))xi{xt^i+,-i) + O(x') ■ (65) 
From this we derive the following linear response 

T , T k 



p{A) = i Yl ^(^*+^)) = Pi"^^ + ^ E E Xi{f-\x,.,))DiA o f-^)if+\£,_,)) + 0{x') (66) 
= p(A)+ / pidx)J2xiix)DiAon{f+\x)) + Oix')- (67) 



Let's consider the case where the memory dependence is given by a sum of memory terms at different times 
in the past: 

oo 

Xt+1 = f{xt) + 'Yxi{xt-i) . (68) 
The response can be derived from Eq. |67|by linearity: 



i = i>0 

In the continuous case, for an equation of the form 



/oo 
p{dx) J2 E Xiix)D{A o f)(f+' (x)) + 0(x'[ 



(69) 



dt 

this response becomes 



d f'°° 

-Xt^Fixt) + xAxt-r) (70) 



6p' 



'T 



/poo poc 
p{dx) ds drxA^)dr(.Mon{rix)) + 0{x'). (71) 
^0 ^0 
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Figure 1: First order response: diagram of the term describing the impact of the averages of the Y variables 
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Figure 2: Second order response: diagram of the term describing the impact of the fluctuations of the Y variables 
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Figure 3: Second order response: diagram of the term describing the memory effect of the X variables on themselves 
mediated by the correlations of the Y variables G\\^ *,6('''ii ''"2)- 
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